In this example we use a very powerful feature of Py-ART: The ability to treat radar gates as a cloud of points and mapping these clouds to a grid. Py-ART uses inverse distance weighted nearest neighbour objective analysis [1] to estimate radar parameters on a cartesion (or other) grid. Py-ART allows arbitary formulations for the sphere of influence. "Under the hood" is a very powerful and flexible technique the uses a KD-Tree to do the searching. While the top level method is for regular rectilinear grids the tree mapper can be used to map to other geometries such as aircraft tracks (in the works). First lets import what we need, load and examine some NEXRAD volumes from the Midwest from a few weeks ago:
In [1]:
import pyart
from matplotlib import pyplot as plt
import numpy as np
from netCDF4 import num2date
from copy import deepcopy
%matplotlib inline
print pyart.__version__
Read in the data from three NEXRAD volumes. To save memory do not read in some of the moments, and only keep the 4 lowest sweeps from each radar. Data available here: Click!
In [2]:
excluded_fields = ['velocity', 'spectrum_width', 'differential_phase']
keep_nsweeps = 4
joliet_radar = pyart.io.read('data/KLOT_110622.nexrad',
exclude_fields=excluded_fields)
joliet_radar = joliet_radar.extract_sweeps(range(keep_nsweeps))
lincoln_radar = pyart.io.read('data/KILX_110437.nexrad',
exclude_fields=excluded_fields)
lincoln_radar = lincoln_radar.extract_sweeps(range(keep_nsweeps))
davenport_radar = pyart.io.read('data/KDVN_110608.nexrad',
exclude_fields=excluded_fields)
davenport_radar = davenport_radar.extract_sweeps(range(keep_nsweeps))
Lets see what fields NEXRAD gives us:
In [3]:
for radar in [joliet_radar, lincoln_radar, davenport_radar]:
print radar.fields.keys()
Lets now plot PPIs of reflectivity ($Z_e$), differential reflectivity ($Z_{dr}$) and correlation co-efficent ($\rho_{hv}$) for the four radars
In [4]:
#set the domain:
max_lat = 43.
min_lat = 40
min_lon = -90.5
max_lon = -86
# create an instance of the display class
display = pyart.graph.RadarMapDisplay(joliet_radar)
# generate the figure
plt.figure(figsize = [17,4])
# first subplot: Reflectivity
plt.subplot(1, 3, 1)
display.plot_ppi_map('reflectivity',
max_lat=max_lat, min_lat=min_lat,
min_lon=min_lon, max_lon=max_lon,
vmin=-8, vmax=64,
lat_lines=np.arange(min_lat,max_lat,.5),
lon_lines=np.arange(min_lon, max_lon, 1),
resolution='i')
# second subplot: differential reflectivity
plt.subplot(1, 3, 2)
display.plot_ppi_map('differential_reflectivity',
max_lat=max_lat, min_lat=min_lat,
min_lon=min_lon, max_lon= max_lon,
vmin=-2, vmax=6,
lat_lines=np.arange(min_lat,max_lat,.5),
lon_lines=np.arange(min_lon, max_lon, 1),
resolution='i')
#third subplot: correlation coef
plt.subplot(1, 3, 3)
display.plot_ppi_map('cross_correlation_ratio',
max_lat=max_lat, min_lat=min_lat,
min_lon=min_lon, max_lon=max_lon,
vmin=.5, vmax=1.,
lat_lines=np.arange(min_lat,max_lat,.5),
lon_lines=np.arange(min_lon, max_lon, 1),
resolution='i')
del display
In [5]:
# create an instance of the display class
display = pyart.graph.RadarMapDisplay(lincoln_radar)
# generate the figure
plt.figure(figsize=[17,4])
# first subplot: Reflectivity
plt.subplot(1, 3, 1)
display.plot_ppi_map('reflectivity',
max_lat=max_lat, min_lat=min_lat,
min_lon=min_lon, max_lon=max_lon,
vmin=-8, vmax=64,
lat_lines=np.arange(min_lat,max_lat,.5),
lon_lines=np.arange(min_lon, max_lon, 1),
resolution='i')
# second subplot: differential reflectivity
plt.subplot(1, 3, 2)
display.plot_ppi_map('differential_reflectivity',
max_lat= max_lat, min_lat=min_lat,
min_lon= min_lon, max_lon= max_lon,
vmin=-2, vmax=6,
lat_lines=np.arange(min_lat,max_lat,.5),
lon_lines=np.arange(min_lon, max_lon, 1),
resolution='i')
# third subplot: correlation coef
plt.subplot(1, 3, 3)
display.plot_ppi_map('cross_correlation_ratio',
max_lat=max_lat, min_lat=min_lat,
min_lon=min_lon, max_lon= max_lon,
vmin=.5, vmax=1.,
lat_lines=np.arange(min_lat,max_lat,.5),
lon_lines=np.arange(min_lon, max_lon, 1),
resolution='i')
del display
In [6]:
# create an instance of the display class
display = pyart.graph.RadarMapDisplay(davenport_radar)
# generate the figure
plt.figure(figsize=[17,4])
# first subplot: Reflectivity
plt.subplot(1, 3, 1)
display.plot_ppi_map('reflectivity',
max_lat=max_lat, min_lat=min_lat,
min_lon=min_lon, max_lon= max_lon,
vmin=-8, vmax=64,
lat_lines=np.arange(min_lat,max_lat,.5),
lon_lines=np.arange(min_lon, max_lon, 1),
resolution='i')
# second subplot: differential reflectivity
plt.subplot(1, 3, 2)
display.plot_ppi_map('differential_reflectivity',
max_lat=max_lat, min_lat=min_lat,
min_lon=min_lon, max_lon= max_lon,
vmin=-2, vmax=6,
lat_lines=np.arange(min_lat,max_lat,.5),
lon_lines=np.arange(min_lon, max_lon, 1),
resolution='i')
# third subplot: correlation coef
plt.subplot(1, 3, 3)
display.plot_ppi_map('cross_correlation_ratio',
max_lat=max_lat, min_lat=min_lat,
min_lon=min_lon, max_lon=max_lon,
vmin=.5, vmax=1.,
lat_lines=np.arange(min_lat,max_lat,.5),
lon_lines=np.arange(min_lon, max_lon, 1),
resolution='i')
del display
Lets to some cleaning up of the ZDR first..
In [7]:
def smooth_diff_refl(radar):
smooth_zdr = radar.fields['differential_reflectivity']['data'].copy()
smooth_zdr[np.where(radar.fields['cross_correlation_ratio']['data']<0.5)] = np.ma.masked
for i in range(smooth_zdr.shape[0]):
smooth_zdr[i,:] = pyart.correct.phase_proc.smooth_and_trim(
smooth_zdr[i,:], 8)
radar.add_field_like('differential_reflectivity',
'differential_reflectivity_smooth',
smooth_zdr)
smooth_diff_refl(joliet_radar)
smooth_diff_refl(davenport_radar)
smooth_diff_refl(lincoln_radar)
Now to do the mapping to a cartesian grid. The pyart.map.grid_from_radars method is well documented here Click!.
In [8]:
mesh_mapped_ALL = pyart.map.grid_from_radars(
(joliet_radar, davenport_radar, lincoln_radar),
grid_shape=(1, 301, 301),
grid_limits=((1000, 1000),
(-300.*1000., 300.*1000.),
(-300.*1000., 300.*1000.)),
grid_origin = (41.5, -88.5),
fields=['reflectivity',
'differential_reflectivity',
'cross_correlation_ratio',
'differential_reflectivity_smooth'],
refl_field='reflectivity',
max_refl=100.,
copy_field_data=False)
lets do some simple investigation of the object returned, a grid object...
In [9]:
print mesh_mapped_ALL.fields.keys()
In [10]:
print mesh_mapped_ALL.axes.keys()
In [11]:
plt.figure(figsize = [15,10])
plt.pcolormesh(
mesh_mapped_ALL.axes['x_disp']['data'],
mesh_mapped_ALL.axes['y_disp']['data'],
mesh_mapped_ALL.fields['reflectivity']['data'][0,:,:],
vmin=-8, vmax=64)
plt.colorbar()
Out[11]:
The ROI field contains the radius of influence used.. in this case Py-ART tried to optimize it for you!
In [12]:
plt.figure(figsize = [15,10])
plt.pcolormesh(
mesh_mapped_ALL.axes['x_disp']['data'],
mesh_mapped_ALL.axes['y_disp']['data'],
mesh_mapped_ALL.fields['ROI']['data'][0,:,:])
plt.colorbar()
Out[12]:
Of course Py-ART contains some nice methods to visualize the data:
In [13]:
max_lat = 43.
min_lat = 40
min_lon = -90.5
max_lon = -86
plt.figure(figsize = [15,10])
display = pyart.graph.GridMapDisplay(mesh_mapped_ALL)
display.plot_basemap(lat_lines=np.arange(min_lat,max_lat,.5),
lon_lines=np.arange(min_lon, max_lon, 1),
resolution='i')
display.plot_grid('reflectivity', vmin=-8, vmax=64)
display.plot_colorbar()
del display
joliet_radar = pyart.io.read('data/KLOT_110622.nexrad')
lincoln_radar = pyart.io.read('data/KILX_110437.nexrad')
davenport_radar = pyart.io.read('data/KDVN_110608.nexrad')
smooth_diff_refl(joliet_radar)
smooth_diff_refl(davenport_radar)
smooth_diff_refl(lincoln_radar)
mesh_mapped_ALL_3d = pyart.map.grid_from_radars(
(joliet_radar,davenport_radar,lincoln_radar),
grid_shape=(35, 301, 301),
grid_limits=((0, 17000),
(-300.*1000., 300.*1000.),
(-300.*1000., 300.*1000.)),
grid_origin=(41.5, -88.5),
fields=['reflectivity', 'differential_reflectivity',
'cross_correlation_ratio',
'differential_reflectivity_smooth'],
refl_field='reflectivity',
max_refl=100.,
copy_field_data=False)
pyart.io.write_grid('data/IL_grid.nc', mesh_mapped_ALL_3d)
In [14]:
mesh_mapped_ALL_3d = pyart.io.read_grid('data/IL_grid.nc')
So now we have a three dimensional data set we can slice and dice to get an idea of the precipitation morphology!
In [15]:
display = pyart.graph.GridMapDisplay(mesh_mapped_ALL_3d)
fig = plt.figure(figsize=[15, 8])
map_panel_axes = [0.05, 0.05, .4, .80]
x_cut_panel_axes = [0.55, 0.10, .4, .30]
y_cut_panel_axes = [0.55, 0.50, .4, .30]
colorbar_panel_axes = [0.05, 0.90, .4, .03]
# parameters
level = 2
vmin = -8
vmax = 64
lat = 41.7092
lon = -87.9820
# panel 1, basemap and radar reflectivity
ax1 = fig.add_axes(map_panel_axes)
display.plot_basemap(resolution='i')
display.plot_grid('reflectivity', level=level, vmin=vmin, vmax=vmax, cmap=pyart.graph.cm.NWSRef)
display.plot_crosshairs(lon=lon, lat=lat)
cbax = fig.add_axes(colorbar_panel_axes)
display.plot_colorbar(label='Eq. Reflectivity Factor (dBZ)', cax = cbax )
# panel 2, longitude slice.
ax2 = fig.add_axes(x_cut_panel_axes)
display.plot_longitude_slice('reflectivity', lon=lon, lat=lat,
cmap=pyart.graph.cm.NWSRef, vmin=vmin, vmax=vmax)
ax2.set_xlabel('Distance from Argonne (km)')
# panel 3, latitude slice
ax3 = fig.add_axes(y_cut_panel_axes)
display.plot_latitude_slice('reflectivity', lon=lon, lat=lat,
cmap=pyart.graph.cm.NWSRef, vmin=vmin, vmax=vmax)
# add a title
slc_height = mesh_mapped_ALL_3d.axes['z_disp']['data'][level]
dts = num2date(mesh_mapped_ALL_3d.axes['time']['data'], mesh_mapped_ALL_3d.axes['time']['units'])
datestr = dts[0].strftime('%H:%M UTC on %d %b %Y')
title = 'Sliced at ' + str(slc_height) + ' meters agl at ' + datestr
fig.text(0.5, 0.9, title)
Out[15]:
In [16]:
fig = plt.figure(figsize=[15, 8])
# parameters
vmin = -1
vmax = 6
# panel 1, basemap and radar reflectivity
ax1 = fig.add_axes(map_panel_axes)
display.plot_basemap(resolution='i')
display.plot_grid('differential_reflectivity', level=level, vmin=vmin, vmax=vmax)
display.plot_crosshairs(lon=lon, lat=lat)
cbax = fig.add_axes(colorbar_panel_axes)
display.plot_colorbar(label='Eq. Reflectivity Factor (dBZ)', cax = cbax )
# panel 2, longitude slice.
ax2 = fig.add_axes(x_cut_panel_axes)
display.plot_longitude_slice('differential_reflectivity', lon=lon, lat=lat,
vmin=vmin, vmax=vmax)
ax2.set_xlabel('Distance from Argonne (km)')
# panel 3, latitude slice
ax3 = fig.add_axes(y_cut_panel_axes)
display.plot_latitude_slice('differential_reflectivity', lon=lon, lat=lat,
vmin=vmin, vmax=vmax)
# add a title
slc_height = mesh_mapped_ALL_3d.axes['z_disp']['data'][level]
dts = num2date(mesh_mapped_ALL_3d.axes['time']['data'], mesh_mapped_ALL_3d.axes['time']['units'])
datestr = dts[0].strftime('%H:%M UTC on %d %b %Y')
title = 'Sliced at ' + str(slc_height) + ' meters agl at ' + datestr
fig.text(0.5, 0.9, title)
Out[16]:
In [17]:
fig = plt.figure(figsize=[15, 8])
# parameters
vmin = -1
vmax = 6
# panel 1, basemap and radar reflectivity
ax1 = fig.add_axes(map_panel_axes)
display.plot_basemap(resolution='i')
display.plot_grid('differential_reflectivity_smooth', level=level, vmin=vmin, vmax=vmax)
display.plot_crosshairs(lon=lon, lat=lat)
cbax = fig.add_axes(colorbar_panel_axes)
display.plot_colorbar(label='Eq. Reflectivity Factor (dBZ)', cax = cbax )
# panel 2, longitude slice.
ax2 = fig.add_axes(x_cut_panel_axes)
display.plot_longitude_slice('differential_reflectivity_smooth', lon=lon, lat=lat,
vmin=vmin, vmax=vmax)
ax2.set_xlabel('Distance from Argonne (km)')
# panel 3, latitude slice
ax3 = fig.add_axes(y_cut_panel_axes)
display.plot_latitude_slice('differential_reflectivity_smooth', lon=lon, lat=lat,
vmin=vmin, vmax=vmax)
# add a title
slc_height = mesh_mapped_ALL_3d.axes['z_disp']['data'][level]
dts = num2date(mesh_mapped_ALL_3d.axes['time']['data'], mesh_mapped_ALL_3d.axes['time']['units'])
datestr = dts[0].strftime('%H:%M UTC on %d %b %Y')
title = 'Sliced at ' + str(slc_height) + ' meters agl at ' + datestr
fig.text(0.5, 0.9, title)
Out[17]: